Energy dynamics in a simulation of LAPD turbulence 

B. Friedman, 1 '!^ T. A. Carter, 1 M.V. Umansky, 2 D. Schaffner, 1 and B. Dudson 3 

^Department of Physics and Astronomy, University of California, Los Angeles, 
California 90095-1547, USA 

^Lawrence Livermore National Laboratory, Livermore, California 94550, 
USA 

^Department of Physics, University of York, Heslington, York YO10 5DD, 
United Kingdom 

Energy dynamics calculations in a 3D fluid simulation of drift wave turbulence in the 
linear Large Plasma Device (LAPD) [W. Gekelman et al, Rev. Sci. Inst. 62, 2875 
(1991)] illuminate processes that drive and dissipate the turbulence. These calcula- 
tions reveal that a nonlinear instability dominates the injection of energy into the 
turbulence by overtaking the linear drift wave instability that dominates when fluc- 
tuations about the equilibrium are small. The nonlinear instability drives flute-like 
(k\\ = 0) density fluctuations using free energy from the background density gradient. 
Through nonlinear axial wavenumber transfer to ku ^ fluctuations, the nonlin- 
ear instability accesses the adiabatic response, which provides the requisite energy 
transfer channel from density to potential fluctuations as well as the phase shift that 
causes instability. The turbulence characteristics in the simulations agree remark- 
ably well with experiment. When the nonlinear instability is artificially removed 
from the system through suppressing k\\ = modes, the turbulence develops a coher- 
ent frequency spectrum which is inconsistent with experimental data. This indicates 
the importance of the nonlinear instability in producing experimentally consistent 
turbulence. 
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I. INTRODUCTION 



It is common practice to study a system's linear stability properties to gain insight into 
turbulent dynamics. It is often easier to calculate and analyze linear modes and growth rates 
than to simulate and analyze nonlinear turbulence. However, there are several situations in 
which linear properties can be misleading in understanding turbulent systems. First, linear 
studies of magnetically confined plasmas that neglect stable branches of the linear disper- 
sion relation often miss details of nonlinear dynamics. For example, stable eigenmodes can 
often impact nonlinear dynamics by providing energy sinks and sometimes energy sources 
not found on the most unstable linear branch^- 10 . Stable eigenmodes can shift the energy 
injection and dissipation ranges, making the turbulent dynamics very different from the Kol- 
mogorov picture of hydrodynamic turbulence^. Second, systems with non-normal modes 
(non-orthogonal eigenvectors) display properties that are unexpected from linear calcula- 
tions—. In fact, systems with non-normal modes even make it difficult to predict dynamics 
when stable eigenmode branches are included in analyses^. Third, linear stability analysis 
can miss crucial nonlinear instability effects, which come in several varieties. 

The most obvious variety of a nonlinear instability effect is that of subcritical turbulence 
in which no linear instabilities exist but turbulence is self-sustained given finite-amplitude 
seed perturbations. Subcritical turbulence is common in hydrodynamics^. While not as 
well-known in plasma physics, several cases of subcritical plasma instabilities have been 
shown in the literature^— . The second variety of nonlinear instability includes cases in 
which a particular linear instability is present in a system, but the turbulence is maintained 
by a nonlinear instability mechanism with different physical origin than the linear instability 
mechanism. This has been explored in tokamak edge simulations in which linear ballooning 
instability drive is overtaken in the saturated phase by a nonlinear drift-wave drive^— . 
Finally, it is often found that a particular linear instability is enhanced, depressed, and/or 
modified in the saturated phase by a nonlinear instability with a similar mechanism as the 
linear instability. In some of these cases nonlinear wavenumber transfers can increase or 
cause drive^^., while in other cases zonal flow effects decrease drive^ 1 ^. 

In order to avoid the pitfalls of relying too heavily on linear stability calculations in form- 
ing conclusions on turbulence characteristics, it is useful to perform turbulent simulations 
and diagnose them with energy dynamics analyses. Energy dynamics analyses track energy 
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input into turbulent fluctuations and energy dissipation out of them. They also track conser- 
vative energy transfer between different energy types (e.g. from potential to kinetic energy) 
and between different scales, waves, or eigenmodes of a system. In all, energy dynamics 
analysis can be used as a post-processing tool to characterize simulation turbulence in order 
to gain insight into underlying physical processes. 

In this study, a simulation of a two-fluid Braginskii model of turbulence in the Large 
Plasma Device (LAPD) is subjected to such an energy dynamics analysis. This reveals that 
a nonlinear instability drives and maintains the turbulence in the steady state saturated 
phase of the simulation. While a linear resistive drift wave instability resides in the system, 
the nonlinear drift wave instability dominates when the mean fluctuation amplitude is over a 
few percent of the equilibrium value. The primary linear instability is the resistive drift wave 
which has a positive linear growth rate for low but finite kn . However, the saturated state of 
the simulated turbulence is strongly dominated by flute-like (k\\ = 0) fluctuations in density 
and potential. The flute-like fluctuation spectrum is generated by a nonlinear instability. 
The nonlinear instability is identified by its energy growth rate spectrum, which varies 
significantly from the linear growth rate spectrum. If k\\ = fluctuations are removed from 
the simulation (while retaining zonal flows), the saturated turbulent state is qualitatively 
and quantitatively different and much less consistent with experimental measurement. 

II. THE DRIFT WAVE MODEL 

A Braginskii-based fluid model 30 is used to simulate drift wave turbulence in LAPD 
using the BOUT++ codo 3 ^. The evolved variables in the model are the plasma density, 
N, the electron fluid parallel velocity v\\ e , the potential vorticity w = Vj_ • (N Q V±4>), and 
the electron temperature T e . The ions are assumed cold in the model (7$ = 0), which 
eliminates ion temperature gradient drive, and sound wave effects are neglected. Details of 
the simulation code and derivations of the model may be found in previous LAPD verification 
and validation studies^"—, although electron temperature fluctuations were not included in 
those studies. 

The equations are developed with Bohm normalizations: lengths are normalized to the 
ion sound gyroradius p s , times to the ion cyclotron time u^ 1 , velocities to the sound speed 
c s , densities to the equilibrium peak density, and electron temperatures and potentials to the 
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equilibrium peak electron temperature. These normalizations are constants (not functions 
of radius) and are calculated from these reference values: the magnetic field is 1 kG, the ion 
unit mass is 4, the peak density is 2.86 x 10 12 cm -3 , and the peak electron temperature is 
6 eV. The equations are: 



d t N = -v E ■ ViV - iVW| n | e + fi N VlN + S N + {0, N}, (1) 

d t v lle = -— ^V„JV - 1.71^V||T e + ^V|,0 - v e v h + {0, % }, (2) 

d t zu = -iVoVpHe - v in w + /i^Vitu + {0, w}, (3) 

d t T e = ~V E ■ VT e0 - lJl^TeoVu^e + T^r^HeVjTe 

lme -v e T e + /i T ViT e + S T + {0, T e }. (4) 



rrii 



In these equations, /iat, ^t, and ji^ are artificial diffusion and viscosity coefficients used 
for subgrid dissipation. They are large enough to allow saturation and grid convergence^, 
but small enough to allow for turbulence to develop. In the simulations, they are all given 
the same value of 1.25 x 10~ 3 in Bohm-normalized units. This is the only free parameter 
in the simulations. All other parameters such as the electron collisionality z/ e , ion- neutral 
collisionality Ui n , parallel electron thermal conductivity K\\ e , and mass ratio ^ are calculated 
from the experimental parameters. There are two sources of free energy: the density gradient 
due to the equilibrium density profile iV , and the equilibrium electron temperature gradient 
in T e0 , both of which are taken from experimental fits. iVo and T e0 are functions of only the 
radial cylindrical coordinate r, and they are shown in Fig. [TJ 

The terms in Poisson brackets are the E x B advective nonlinearities, which are the 
only nonlinearities used in the simulations. The numerical simulations are fully spatial in all 
three dimensions (as opposed to spectral) and use cylindrical annular geometry (12 < r < 40 
cm). The radial extent used in the simulation encompasses the region where fluctuations are 
above a few percent in the experiment. The simulations use periodic boundary conditions 
in the axial (z) direction and Dirichlet boundary conditions in the radial (r) direction for 
the fluctuating quantities. 

Simulations also use density and temperature sources (S n and St) in order to keep the 
equilibrium profiles from relaxing away from their experimental shapes. These sources 
suppress the azimuthal averages (m = component of the density and temperature flue- 



tuations) at each time step. The azimuthal average of the potential <fi is allowed to evolve 
in the simulation, allowing zonal flows to arise. The parallel current, which is often found 
explicitly in these equations is replaced here by J\\ = —N v\\ e . 
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FIG. 1. The profiles of density Nq and electron temperature T c q used in the simulations normalized 
to their peak values of 2.86 x 10 12 cm -3 and 6 eV, respectively. 

Some basic statistical properties of the density fluctuations of the simulation are shown 
in Fig. [2] and are compared to the corresponding results from the experiment on which 
this simulation is based. The simulation reproduces these characteristics of experimental 
measurements with rather good qualitative and quantitative accuracy. 

III. ENERGETICS MACHINERY 

In order to perform an energy dynamics analysis on the simulation, expressions for the 
energy and energy evolution must be derived from Eqs. [T]|H To start, an expression for the 
normalized energy of the wave fluctuations in the drift wave model is defined as: 

E=U{N* + \Tl + ^vj e + N (V ± cf>r)dV. (5) 
i jy a rrii 

The N 2 contribution is the potential energy due to density fluctuations, |T e 2 is the electron 
temperature fluctuation potential energy, ^v 2 e is the parallel electron kinetic energy, and 
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FIG. 2. a) The power spectral density of the density fluctuations, showing the results from simu- 
lation versus experiment, b) the probability distribution function of the density fluctuations, and 
c) the RMS amplitude of the density fluctuations as a function of radius. 



iVo(V_i_0) 2 is the E x B perpendicular kinetic energy. These energy-like expressions are 
defined in this way so that they are conserved individually by their respective advective 
nonlinearities. While the physical energy contains extra factors of N Q and T e0 , the physical 
energy does not preserve the property of conservative nonlinearities of Eqs. HHand therefore 
produces a more complicated analysis. Analyzing an energy-like expression such as that in 
Eq. [5l however, can be just as illuminating, and in this case simpler, than analyzing the 
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physical energy. The expression in Eq. [5] will henceforth be referred to simply as the energy. 

Furthermore, it is often most instructive to analyze the spectrally decomposed energy 
dynamics. To do this, each fluid field (N, T e , v\\ e , 0) at a given time is Fourier decomposed as 
F(r,8,z) = Xlfc /fe( r ) e ^ m6 ' +fczZ ' ) ) where the subscript k represents the spectral wavenumbers, 
(m,n). m is the azimuthal wavenumber while n is the axial integer wavenumber such that 
k z = k\\ = 2ixn/l z . Note that the radial direction is not spectrally decomposed because the 
radial dependence of the profiles and differential operators complicates the analysis. With 
this, the energy of each Fourier k = (m, n) mode is 



Etot(k) = I ( kh 2 + 1|^| 2 + ^KH 2 + No 



Or 



+ ^|<tf>, (6) 



where the brackets () represent the radial integral: J^ b rdr. The energy evolution for each 
Fourier mode of each field has the form: 

= Q 3 {k) + C,{k) + D,(k) + £ k'). (7) 

%' 

The index j stands for each field, (n,t,v,(p), and the sum over j gives the total energy 
evolution. Note that with the conventions used, the symbol n denotes both the axial mode 
number as well as the Fourier coefficient of the density fluctuation. The differences should 
be clear in context. The derivation of Eq. [7] is given in the Appendix along with the full 
expressions for each of the parts. Tj(k,k') is the nonlinear energy transfer function that 
comes from the advective nonlinearities. It describes the nonlinear energy transfer rate of 
modes k' = (W, n') and k — k! = (m — m', n — n') to the mode k = (m, n). In other words, 
a positive value of Tj(k,k') indicates that fluctuations at wavenumber k gain energy from 
gradient fluctuations at wavenumber k! and flow fluctuations at wavenumber k — k! . When 
summed over k! as in Eq. [?l the result is the total nonlinear energy transfer into mode k. 
Note that ^likk' Tj(k> k') = because the nonlinearities conserve energy individually in each 
of Eqs. [10 This is easily proven by the following identity: 

/ q{p,q}dn= [ p{p,q}dn = 0, (8) 
Jn Jn 

which holds when boundary conditions are periodic or zero value as they are in the 
simulation. The fact that the advective nonlinearities conserve energy means that they 



can transfer energy between different Fourier modes but they cannot change the energy of 
the volume-averaged fluctuations as a whole. Only the linear terms can change the total 
energy of the fluctuations. Other possible nonlinearities that do not conserve energy are 
not included in the model equation set or in the simulations for simplicity of the energy 
analysis. Furthermore, it is convenient for the simulations to employ an energy conserving 
finite difference scheme for the advective nonlinearities to reflect this analytic property of 
the equations. However, most common numerical advection schemes do not conserve energy 
for finite grid spacing. Therefore, an Arakawa advection scheme 36 that conserves energy of 
the advected quantity is used for the nonlinear advection terms in the simulations. 

The linear terms in Eqs. dHdo not conserve energy individually or as a whole. The linear 
terms are broken up into three contributions in Eq. [7J Dj(k) represents nonconservative 
energy dissipation due to collisions, artificial diffusion and viscosity, and the density and 
temperature sources. Each contribution to Dj{k) is negative, and the exact expressions are 
given in the Appendix. Cj(k) contains the linear terms dubbed "transfer channels"—. They 
are rewritten here: 



C n (k) = Re {(-ik z N vpil) } (9) 

C v {k) = Re {{-ik z N np)l + ik z N cf>pjt - 1.71^ 2 T e0 t^)} (10) 

C4(k) = Re{(ik t N v&l>l)} (11) 

C t {%) = Re{(-1.7Uk z T e0 v,:tl)} (12) 

Notice that C n (k) + C v {k) + C${k) + C t {k) = 0, which is most clearly seen upon conju- 
gation of C v (k) inside the real part operator. This is the reason why these terms are called 
transfer channels. They represent the transfer between the different types of energy of the 
different fields (N,<p,T e -h- v\\ e ), but taken together, they do not create or dissipate total 
energy from the system. The only energy field transfer in this system occurs through the 
parallel electron velocity (parallel current) dynamics. There is no direct transfer between 
the state variables N,<p, and T e . Altogether, the coupling through the parallel current is 
called the adiabatic response. It is an essential part of both the linear and nonlinear drift 
wave mechanisms^ 3 - 1 ^. The adiabatic response moves energy from the pressure fluctuations 
to the perpendicular flow through the parallel current. 
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Finally, the Qj(k) terms represent the nonconservative energy sources. They are rewritten 
here: 



— * / / 777 

Q n (k) = Re\( drNo^nt 



Q v (k) = Re < ( ik z 




■npjt + ik z (l - N )4>pjl + 1.71ik z (T e0 - l)tp;l 



(13) 



qS) = o 




(15) 
(16) 



Q n {k) is the energy extraction from the equilibrium density profile into the density fluc- 
tuations. This term may have either sign depending on the phase relation between <f)^ and 
nj:, so it can in fact dissipate fluctuation potential energy from the system as well as cre- 
ate it at each k. Qt{k) is completely analogous to Q n (k) but for the temperature rather 



these terms are the equilibrium gradients, which is evident because if the profiles were flat 



expression analyzed here Q v (k) is small compared to the parallel velocity dissipation D v (k), 
meaning that there is no net energy entering the system directly through the parallel veloc- 
ity as expected. The energy expression used here is therefore a good proxy for the physical 
energy with the added benefit that the expression used here conserves the nonlinearities. 

IV. NONLINEAR ENERGY DYNAMICS 

A. Energy Spectra 

Figure [3] shows the time evolution of the total energy of the fluctuations. The simulation 
starts with a random initial perturbation, and the fluctuations grow exponentially due to 
the linear drift wave instability until the energy level reaches about 0.01, where the energy 
is fairly equally divided between n = and n = ±1 modes. Then, the nonlinear instabil- 
ity takes over and the fluctuation energy continues to grow until reaching saturation. All 



than the density. Q v (k) is parallel kinetic energy extraction or dissipation. The sources of 




9 



0.30 
0.25 
0.20 
0.15 
0.10 
0.05 
0.00 








turbulent stage 

( E tot)V 



100 150 
Time Step 



200 250 



FIG. 3. Time evolution of the volume-averaged total energy. Each time step is 400/w c j ~ 170fis 
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FIG. 4. a) E n (k), b) Et(k), c) E^(k), and d) E v (k) in the m — n plane averaged over time during 
the saturated turbulent phase. Note the different scales used on each figure and also that the 
vertical white lines at m = in a) and b) are due to the density and temperature sources which 
subtract out this component of the fluctuations. The zonal flow, defined as the n = 0, m = 
component of the E x B velocity, has an energy of magnitude of 1.4 x 10 -4 as seen in c). 
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analysis shown below is done by time averaging over the saturated (turbulent) stage. The 
turbulent spectral energy, defined in Eq. [6] is shown in Fig. HJ The energy is broken up into 
its different types (e.g. perpendicular kinetic energy: E^). There are a few clear nonlinear 
properties seen in these figures. The first is that the energy is located in different spectral 
regions for the different energy types. This has to be a nonlinear effect because the linear 
eigenmodes are Fourier modes in the azimuthal and axial directions and all fields grow at 
the same rate for an eigenmode. Another property unexpected from linear stability analysis 
is that most of the potential and perpendicular kinetic energy (E n , E t , and E^) is contained 
in n = (/c|| = 0) structures, which are often called flute modes. Previous studies pointed 
out this flute mode dominance in LAPD simulations^ 1 ^. The study by Rogers et al.— , how- 
ever, used a momentum source that produced a large radial electric field, possibly leading 
to a dominating Kelvin-Helmoltz instability at k\\ = 0. Such a feature is unexpected in this 
study because there is no n = linear instability present in the system, which is confirmed 
by eigenvalue calculations^. The only linear instability of the system is the linear resistive 
drift wave instability, which requires finite n to provide the phase shift and state variable 
coupling to drive the waves unstable. Perhaps equally unexpected is the complete lack of 
parallel kinetic energy in the n = range. The E v spectrum looks like a traditional linear 
drift wave spectrum, but does not match the other fields, which is atypical of linear drift 
waves. 

B. Description and Evidence for the Nonlinear Instability 

The flute mode dominance has to be a nonlinear effect because linear drift waves re- 
quire finite n. However, the cause of the flute dominance is not simple cascade dynamics, 
a secondary instability, nor a flow-driven flute-like instability such as Kelvin-Helmholtz or 
interchange. Rather, the cause is a primary nonlinear instability as has been reported in 
previous simulations of plasma edge turbulence^^. This nonlinear instability is a multi-step 
process that is outlined in Fig. EJ In the first step, n = density and potential fluctua- 
tions nonconservatively draw energy from the equilibrium density gradient as prescribed by 
Q n {m, n = 0) defined in Eq. [13], and feed this energy into the n = density fluctuations only. 
The nonconservative linear terms, after all, can only inject, dissipate, or transfer energy at 
one wavenumber at a time, so it takes n = fluctuations to nonconservatively inject energy 
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FIG. 5. Diagram of the nonlinear instability process that drives flute modes. 



into n = fluctuations. Note also that the temperature fluctuations work in the same way 
as the density fluctuations, and one could replace n(m, n) (the spectral density component) 
with t(m,n) (the spectral temperature component) in the diagram. The density and tem- 
perature are analogous and work in parallel, however the temperature fluctuations are a 
few times smaller than the density fluctuations and provide weaker drive because parallel 
heat conduction strongly dissipates the temperature fluctuations. This is why the diagram 
highlights the density contribution (this observation is consistent with other work in edge 
turbulence^). This nonconservative injection does not occur for infinitesimal perturbations; 
a finite-amplitude n = seed perturbation is required. In the simulations, this seed is pro- 
vided by nonlinear transfer from n = 1 drift wave fluctuations which dominate the linear 
phase of the turbulence simulation. 

In the second step of the diagram, these n = density fluctuations conservatively transfer 
energy to n ^ density fluctuations by the nonlinear T n (k, k') transfer process. The third 
step involves the transfer at finite n from the density fluctuations to the potential fluctuations 
by way of the parallel current in the adiabatic response. The Cj(k) terms describe this 
adiabatic response. Fourth and finally, the T^(fc, k') interaction conservatively transfers 
energy from n ^ to n = potential fluctuations in inverse fashion, providing the 
necessary potential flute structures for the first step. 

The evidence for the dominance of this mechanism is best shown with help from the 
energy dynamics machinery derived in section IIHI Figure [6] summarizes the effects of the 
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FIG. 6. a) The solid curves quantify the energy dynamics of the density potential energy averaged 
over time during the saturated turbulent phase. The notation n ± 1 represents the summation over 
the n = 1 and n = — 1 curves, b) The energy dynamics of the temperature potential energy, c) the 
perpendicular kinetic energy, and d) the parallel dynamics (adiabatic response). The contributions 
to C v {k) in d) are broken up with C v (k) = —C n {k) — Ct(k) — CAk). The density and temperature 
sources are neglected in a) and b) respectively. They only contribute at m = 0. 



nonconservative linear terms, which are fully responsible for injecting energy into the fluc- 
tuations. Figure [6fa) shows the E n dynamics separated into different parallel wavenumbers 
and plotted against the azimuthal wavenumber m. Clearly, most of the energy is injected 
into n = density structures, while only a small amount of energy is injected into n — ±1 
structures despite the fact that the linear instability is active only at n ^ 0. The large 
positive Q n + D n (injection plus dissipation) at n = provides evidence for the first step 
of the diagram in Fig. [51 Note however that the dissipation from the source, which acts 
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entirely at m = 0, is neglected in this figure because it is so large (about 5 x 10 -5 ) that 
it would compress the other lines too much. Modes with |n| > 2, on the other hand, play 
a negligible role in density injection, dissipation, and transfer. Furthermore, all of the net 
energy injected into the density fluctuations (Q n + D n ) is transferred out (C n ) to the parallel 
current (electron velocity), which only occurs at finite n, almost entirely at n = ±1. The 
net change of E n , which is the sum Q n + D n + C n over all m and n is approximately zero 
because this analysis is averaged over the steady state turbulence, although this is not so 
evident in Fig. [6^ without the source dissipation. The necessary balance implies, as will be 
proven later, that the nonlinearities transfer energy from n = to n = ±1 modes, where 
that energy can then be transferred to the parallel current. 

Figure E](b) shows the temperature potential energy dynamics. Again flute structures 
inject energy into the fluctuations, but unlike in the density case, n = ±1 modes dissipate 
more energy than they inject. Moreover, the small value of Ct reveals that the temperature 
fluctuations inject only a small amount of energy into the parallel current compared to the 
density fluctuations. Despite the fact that the equilibrium temperature gradient is nearly as 
steep as the density gradient at its steepest point, its free energy is not used efficiently by the 
waves in the sense that it is largely dissipated before being transferred to the electrostatic 
potential. The reason for the difference between the density and temperature responses 
is the extra dissipation routes for the temperature fluctuations, namely, the parallel heat 
conduction and electron- ion heat exchange. One should therefore be careful in assuming 
that adding free energy sources to an analysis will automatically increase instability drive. 
The same type of result was seen in a study of tokamak edge turbulence^, although there, 
the temperature fluctuations were even more dissipative than in this study in that they 
actually drew energy from the parallel current. 

Next, Fig. [6]Jc) illustrates the perpendicular kinetic energy dynamics provided by the 
electrostatic potential 0. Since there is no free energy source for the potential (Q^ = 0), the 
potential fluctuations derive their energy from the parallel current through the transfer 
channel, which is positive everywhere and only finite for finite n. Otherwise, ion-neutral 
collisions and viscosity dissipate energy from the potential fluctuations as shown by the 
curves. An interesting detail seen in this figure is that modes with |n| > 1 actually 
contribute to the transfer channel and dissipation, whereas these modes are negligible for 
the other fields. 
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The last piece, the parallel dynamics, also called the adiabatic response, is displayed in 
Fig. [6(d). The primary effect of the adiabatic response is to take energy from the density 
fluctuations and transfer it to the potential fluctuations, which only occurs at finite parallel 
wavenumber. This effect corresponds to the third step in Fig. Moreover, resistivity dis- 
sipates a substantial portion of this parallel kinetic energy and, although not evident from 
this figure, provides the primary phase shift mechanism between the density and potential 
that allows for instability. The temperature fluctuations also provide energy to the poten- 
tial fluctuations through this response (— C t ), although it is much weaker than the density 
fluctuation route. 
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FIG. 7. a) Conservative nonlinear energy transfer functions Tj(k, k') summed over k' and n. The 
line labeled cj) represents the function n T ( f,(k, k'), which is a function of m. The line labeled 
is the total energy transfer (the sum over the other lines), b) Transfer functions summed over 
k' and m. Note that T v (k,k') is multiplied by 20 in both figures to make it visibly different from 
zero. 

Steps 2 and 4 in Fig. [5] come not from the nonconservative linear terms in the equations, 
but from the conservative nonlinear advective terms. The interactions described by the 
advective nonlinearities are in the nonlinear transfer functions: Tj(k,k r ) in Eq. [7J It is 
difficult to study the Tj(k,k') functions because they are four dimensional functions of 
(m,n,n' ,m'), which makes visualization challenging. It is therefore convenient to sum over 
various transfers or look at specific wavenumber transfers of interest. The most easily 
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decipherable results that complement the results of Fig.[6]are shown in Fig. [7J First, Fig. UJa) 
sums the transfer functions over (n,m',n'), leaving only a function of m, which illustrates 
the aggregate azimuthal mode number transfers. Note that the sum of each individual 
curve over m is zero because the nonlinearities are conservative. Positive values in Fig. [TJa) 
indicate energy transfer into structures with azimuthal mode number m, while negative 
values indicate energy transfer out of structures with corresponding mode number m. The 
density and temperature nonlinearities are qualitatively similar in that they cause both 
forward and inverse transfer out of the wavenumbers that nonconservatively inject the most 
energy. The potential (polarization) and parallel velocity nonlinearities cause inverse transfer 
to low wavenumbers. 

Figure UJb) displays the conservative transfers summed over (m,m',n'), leaving only 
a function of n, which describes transfer into and out of different parallel modes. This 
is the figure which provides evidence for steps 2 and 4 of the instability diagram. Now, as 
expected from step 2 of the diagram and foreshadowed by Fig. E(b), density potential energy 
is aggregately transferred from n = flute modes into n ^ modes, specifically n = ±1. 
This can be called a direct transfer in analogy with the terminology used for hydrodynamic 
wavenumber cascades. The temperature fluctuations have the same transfer trends as the 
density fluctuations, while the parallel velocity exhibits direct transfer, although from \n\ = 1 
to higher modes since there is never any n = energy in the parallel velocity. The potential 
fluctuations, on the other hand, exhibit inverse parallel wavenumber transfer (step 4 of the 
diagram), populating n = potential structures. This nonlinear transfer is the only way to 
drive energy into n = potential structures because = 0. That completes the evidence 
for the nonlinear instability picture along with further details of both the conservative and 
nonconservative energy dynamics. 

C. The Global Energy Injection and Dissipation Picture 

The details of the energy dynamics given above are important but can obscure the most 
significant results. Specifically, Fig. [6] contains a lot of details that can be contracted by 
summing over the different energy types. Figure |8] does this, showing the total spectral 
nonconservative energy dynamics. Figure Eta), which is a plot of the nonconservative rate 
of change of the total energy, 9E ^ \ NC = J2jQj(k) +Cj(k) + Dj(k), reveals a global picture 
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FIG. 8. a) The total spectral nonconservative energy injection ^ - \ NC and b) the spectral 
nonconservative growth rate spectrum 7t(&) of the turbulence compared to the linear growth rate 
spectrum, 7l(/c). The solid lines represent 7t(&) which is calculated using Eq. [T71 averaged over 
the saturated turbulent phase, while the dashed lines represent jL^k) and are calculated with the 
same equation using the linear phase of the simulation. 

in wavenumber space of where the total energy is injected into the system and where it is 
dissipated. Namely, energy injection occurs on average at (n = 0,3 < m < 45), while it is 
dissipated everywhere else including at n = ±1 for all m. It is obvious that the nonlinear 
wavenumber transfers must take energy from the injection region to the dissipation region on 
average, and that is consistent with what was shown in Fig. [71 This further reveals a picture 
quite different than what one would expect from the standard picture of plasma turbulence 
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in which energy is injected where the linear growth rate is positive and dissipated where it 
is negative. The picture here is quite the opposite - energy is injected where there is no 
linear instability and dissipated in part where there is one. To clarify this point, the linear 
growth rate 7x(fc) versus turbulent growth rate 7r(&) spectra are shown in Fig. Efb). The 
growth rates are calculated using: 
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{k)= (^Q^+C^ + D^/^^E^. (17) 

The turbulent growth rate spectrum, simply means that Eq. [T7] is calculated using the 
terms from the saturated turbulent phase of the simulation. Note that the linear growth 
rate is positive for (n = —1,35 < m < 95) and negative everywhere else. The turbulent 
growth rate is positive only for (n = 0, 3 < m < 45). The linear and turbulent spectral 
injection regions do not even overlap. Seemingly, the linear physics is completely washed 
out in the turbulent state. 



V. LINEAR VERSUS NONLINEAR INSTABILITY DRIVE 




FIG. 9. a) The turbulent growth rate spectrum 7t(&) with n = density and potential components 
removed from the simulation compared to the linear growth rate spectrum 7z,(A;). b) The solid 
lines are the same 7t(&) spectrum as the solid lines in (a), but the dashed lines are the turbulent 
growth rate spectrum when the zonal flows are retained in the simulation. 
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Although the nonlinear flute mode dynamics present a clear case of nonlinear instability, 
the ii^0ffn = energy path is not a necessary feature of nonlinear drift wave instabilities, 
which is clear in tokamak studies of drift wave turbulence^ - —. Furthermore, the periodic 
axial boundary conditions used in the LAPD turbulence simulation are obviously unphysical, 
and more realistic boundary conditions may change the parallel dynamics disallowing an 
exact n^0f>n = path. 

In essence, it is interesting to test the robustness of nonlinear instability in this system. 
In particular, how important are the idealized flute modes to the nonlinear instability? 
They are after all, not essential to the otherwise similar nonlinear drift-like instabilities 
in the tokamak edge simulations^ - —. Now, there are a few ways to eliminate the flute 
modes in the simulation such as eliminating one of the nonlinearities that is essential to the 
nonlinear instability process described in Fig. [5j However, simply subtracting out the n = 
components of the density and potential at each simulation time step retains more of the 
physics that may be essential to cause nonlinear instability. The energy dynamics of such 
a simulation, which are succinctly summarized by the growth rate spectrum, are shown in 
Fig. [9j Interestingly, with no n = fluctuations, the turbulent growth rate spectrum 7t(&) 
is nearly identical to the linear growth rate spectrum 7l(/c), as seen in Fig. [9]^a). It is noted 
that subtracting out the n = potential component removes the zonal flow (m = n — 0) 
from the system, providing a possible explanation for the large change in behavior of the 
turbulent growth rate spectrum. However, this hypothesis is dispelled by the analysis of 
a simulation in which only the n = 0, m ^ potential components are removed while the 
zonal flow is left intact. The growth rate spectrum of this simulation, shown in Fig. [9b, 
reveals that the zonal flow plays a minimal role in the nonlinear instability dynamics. The 
zonal flow simply decreases the growth rates by a small amount, causing no change to the 
qualitative picture. 

The flute modes are therefore necessary for a nonlinear instability to overtake the linear 
instability in driving the turbulence, making this qualitatively different from the tokamak 
edge studies. Furthermore, it is clear from this result that a 2D simplification using a fixed 
parallel wavelength like the 2D Hasegawa-Wakatani model^ does not support a nonlinear 
instability. Nevertheless, one indication that the nonlinear flute-driven instability is im- 
portant in reproducing experimentally consistent turbulence is that the turbulence of the 
simulation with the n = components removed becomes overly coherent. This can be seen 

19 



1 n _ 5 
1U 






— Experiment i 


1 rv^ 
1U 


— ^^^^^^^^^^^ C i tyi ii 1 a on - 

_^ 011I1U.1CI L1U11 




** J"^^ — No n=0 


1 rr^ 


^^^^^ 11 U 11 — \J 










10" 8 








10" 9 








10 10 








10" 11 








1Q -12 


'■1 ' ■ ' v' 



10 3 10 4 10 5 

f (Hz) 

FIG. 10. Comparison of the frequency spectra. Notice the spectra with n = components removed 
is not broadband, but has a clear peak, which is inconsistent with experiment. 

in the frequency spectrum, which is shown in Fig. [TU] compared with the experimental spec- 
tra and the spectra of the standard nonlinear-instability-dominated simulation. While the 
standard simulation with the n = modes retained compares well with experimental data, 
the spectrum of the simulation with n = modes removed does not. It is not broadband, 
having a large peak, which is inconsistent with experimental spectrum. Apparently the 
n = modes and the nonlinear instability are important for reproducing experimentally 
relevant turbulence. A more direct test with realistic axial boundary conditions is left for a 
future study. 



VI. CONCLUSION 

In contrast to experiments, simulations provide vast quantities of spatial information, 
and can therefore be useful in illuminating physical processes responsible for driving and 
saturating turbulence. It is possible to get more than fluctuation levels, flux values, diffu- 
sivities, and spectra from simulations. The kind of energy analysis used in this study is one 
way in which detailed physics can be drawn from a turbulence simulation. Here, energy 
dynamics analysis shows a complex picture of turbulent energy injection, transfer, and dissi- 
pation. Such a picture was certainly not evident a priori. Other more advanced procedures 
such as eigenmode decompositions^ or proper orthogonal decompositions (POD) 10 , which 
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are extensions of this procedure, can reveal even more physical processes, especially those 
involving saturation. 

In this study, a partial spectral decomposition and energy dynamics analysis was suf- 
ficient to reveal the dominance of a nonlinear instability in driving and maintaining the 
turbulence. The nonlinear instability works by driving kn = pressure fluctuations using 
k\\ = pressure and potential fluctuations to access the free energy pressure gradients. These 
k\\ = pressure fluctuations are not driven by conservative nonlinear energy transfer from 
linear drift waves nor by some primary linear flow-driven instability. The k\\ = potential 
fluctuations are driven through the finite kn adiabatic response in combination with forward 
and reverse axial wavenumber transfers. Not only does the nonlinear instability require the 
k\\ = fey 7^ transfer path to operate, but the simulation requires this to produce ex- 
perimentally consistent broadband turbulence. In the future, this study will be expanded to 
include different and more realistic axial boundary conditions - including conducting plate 
boundaries - to further test the importance of the flute modes in creating broadband turbu- 
lence. Furthermore, different operating conditions in LAPD, including those that produce 
large mean flows will be simulated and studied to test for the emergence of new dominant 
instabilities, possibly nonlinear ones. 

Understanding nonlinear instabilities is important because they can invalidate the use 
of quasilinear flux estimates and linear mixing length arguments of turbulent transport 
levels when linear instabilities are insignificant in the turbulent state. Simple rules for when 
nonlinear instabilities will act or overtake linear ones are needed, and one attempt at such 
a rule has been made elsewhere for drift wave turbulence^. That rule states that nonlinear 
instabilities will overtake linear instabilities when < which is true for the parameters 
used in this study. However, more study of this rule and others is warranted, and will 
be important as long as linear calculations are used to inform predictions of turbulence. 
Nevertheless, full nonlinear simulations and energy dynamics analyses are most informative 
and should be used to obtain details of plasma turbulence mechanisms. 
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Appendix A: Explicit Calculation of the Energy Evolution of a Fourier Mode 

The energy evolution for each Fourier mode can be obtained by Fourier decomposing 
each of Eqs. [T]-H]and then multiplying the density, electron parallel velocity, vorticity, and 
electron temperature equations by the complex conjugates of the density, velocity, potential, 
and temperature respectively, and integrating over space. Adding the resulting equations 
gives the energy evolution of each Fourier mode. 

Take the density equation as an example for this procedure. The decomposition for the 
density is: 



N(r,9,z,t) =Y,n i :{r,t)e i( - me+k * z \ 



(Al) 



Recall that the subscript k is short for (m, n) as the decomposition is a 2D Fourier 
decomposition in the azimuthal and axial directions, making the sum over k truly a double 
sum over m and n. Furthermore, positive and negative m and n are included in the sums 
to ensure reality of N, which also requires that = n%. Similar decompositions are used 
for v\\ e and <fi. The density source is azimuthally symmetric, so it is decomposed as: 



S N {r,z,t) = ^2S Nkz (r,t)e 



ih.-.y. 



(A2) 



Substituting this decomposition into Eq. CD gives: 



Qnt „■ 



k_ i(m9+k z z) 



dt 



vm 



m 



— d r N (f)j: - ik z N vj: + ii N (d;n% + -d r n^ - 



i(m8+k z z) 



k,k' 



i(m+m')0+i(k z +k' z )z 



(A3) 
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Multiplying through by n*>„e irn " e %k " z and integrating over space (and permuting primes) 
gives: 



dn 



k * 

n, 



dt k 



(SNk z n* m=okz ) 



+ ^ iw!n v d rh-v n \ m')d r n %l (t) % _ %l nl) ^ . (A4) 

Finally, taking the real part of this equation results in: 



ld\nj: 



2 dt 



= Re{(S Nkz n* m =o,k z )} 



+Re | ( r ? ( im ' n k' d ^k-k' n *k ~ ^ m ~ m ') d r n k'h-k' n l) ^ | • ( A5 ) 

Note that taking the real part of the equation produces the expected energy-like term on 
the left hand side because: 

i%!! =JJe /^f„a. ( A6) 

2 dt \ dt k J v ; 

Breaking the result into explicit parts: 



= Q n (k) + C n (k) + D n (k) + Y,Tn(k, k') (A7) 

E n (k) = l(\n n \ 2 ) (A8) 

Q n (k) = Re^-^drNo^n*^ (A9) 

C n (k) = Re{(-ik z N v^nl)} (A10) 

£>„(£) = Re |^/iv(<9 r 2 ng + ^ r n^ - ^-n^n^ + ^fc z < =0 , fc ^| (All) 

T n (fc, fc') = i?e | [im'n^dr^k-k^l ~ K m ~ m ') d r n k^k-k' n l) ^ } ( A12 ) 
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Q n (k) is the energy injection, C n (k) is the transfer channel, D n (k) is dissipation, and 
T n (k, k') is spectral energy transfer. The same type of procedure may be applied to Eqs. [2]- 
IU However, the double primed conjugate multiplications (as in the step between Eqs. IA3I 
and IA4I) must be done with the Fourier fields, ^vp,, —<fip,, and |tg„ rather than vp,, mp, 
and tfc,. These produce the correct energy terms, and most importantly still conserve the 
nonlinearities. The corresponding expressions for the perpendicular kinetic energy are: 



= Q^k) + C^k) + D^k) + P) 



(A13) 



E^k) 



Or 



m 



+ ^o — |^i 



Q4>(k) = 

Ct(k) = Re{(ik z N vtff>l)} 



D^k) = Re M-^^-oj^ + -d r wj: - -^^i)<t>\ - t> in E^{k) 
T^k, k') = Re l(~ ^ w k' d r^~k'^\ ~ *( m _ ^dr^p^k-ft^l) 



(A14) 

(A15) 
(A16) 

(A17) 
(A18) 



and for the electron temperature potential energy: 



-JP = Q t (k) + C t (k) + D t (k) + Tt(k, k') (A19) 

ft 

E t {k) = \{U 2 ) (A20) 

Q t {k) =ite|/-~ft.Teo^jU (A21) 

C t ik) = Re{{-1.71ik z T e0 v^q)} (A22) 

T t (k, k')=Re{/^- {im'tpdr^ptt - i(m - m^drtp^q) ) ) (A24) 



and for the parallel kinetic energy: 
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— ^ = Q v (k) + C v (k) + D v (k) + J2 T ^ fc') (A25) 

k' 

E ^) = \^M 2 ) (A26) 
g w (fc) = ite | (ik j^'^ n^ + ik z (l - N Q )<t> u vl + 1.7HA; 2 (T e0 - 1)^ | (A27) 

C^fc) = Re { (-ih z N ngv* n + ik z N Q <f>p}\ - 1.7H^T e0 ^) } (A28) 

A,(fc) = i?e j^-z/ e ^KH 2 ^} (A29) 

T„(fc, fc') = Re { ^ / - {im'vpdr^pvl - i(m - m')d r v p ^_ p vl) \ \ (A30) 



The transfer channel C v (k) is specifically set so that C n (k) + C t (k) + C^{k) + C„(fc) = 0. 
The source Q v {k) is the left over quantity, which can have any sign and contributes to the 
overall energy evolution. 
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